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1.  Introduction 

Oceanic  and  atmospheric  processes  are  characterized  by 
motions  of  different  scales  which  often  manifest  themselves 
by  distinct  maxima  in  the  spectral  density  of  a  background 
model  state  and  in  the  spectra  of  background  error 
covariance  (BEC)  matrices.  Numerical  modelling  of  the 
multi-scale  BEC  structures  in  variational  data  assimilation 
is  a  challenging  task  which  has  recently  drawn  a  significant 
attention  due  continuous  growth  of  computer  power.  In 
particular,  it  is  desirable  to  formulate  the  scale-dependent 
BEC  operators  (Dance,  2004),  which  can  account  for 
smaller-scale  components  present  in  very  high -resolution 
models.  The  term  ‘multi-scale  correlation  function’  is  also 
used  in  the  theory  of  turbulence  and  reflects  covariances  of 
multifractal  nature  characterized  by  the  power  law  decay  of 
correlations  (e.g.  Mandelbrot,  1997). 

A  straightforward  way  to  construct  multi-scale  BEC 
operators  is  to  use  suitable  superpositions  of  the  single¬ 
scale  correlation  functions  for  modelling  the  BEC  matrix 
elements  (e.g.  Hristopulos,  2003;  Gaspari  etai,  2006).  In 
this  approach,  the  resulting  spectrum  is  difficult  to  control 
directly  by  the  free  parameters  of  the  correlation  functions 
and  care  should  he  taken  to  maintain  positive  definiteness 
of  the  correlation  matrix. 

A  promising  approach  is  to  introduce  scale  separation 
in  the  BEC  models  by  splitting  the  covariance  matrix  into 
several  additive  single-scale  components  (e.g.  Wu  etai , 
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2002;  Purser  et  aL,  2003)  and  perform  assimilation  on  a 
sequence  of  grids  with  increasingly  fine  resolution  (Li  et  aL , 
pers.  comm.  2012). 

A  multi-scale  BEC  operator  can  also  he  constructed 
using  a  polynomial  of  the  discretized  diffusion  operator 
for  representing  the  inverse  covariance.  This  approach  has 
been  studied  by  many  authors  (e.g.  Sasaki,  1970;  Wahba 
and  Wendelberger,  1980;  Purser,  1986;  McIntosh,  1990;  Xu, 
2005).  Its  attractive  features  are  the  flexibility  in  controlling 
the  BEC  spectrum  and  the  low  cost  of  computing  the  action 
of  the  inverse  BEC  matrix  on  a  state  vector.  In  practice, 
however,  applications  of  this  approach  were  limited  to 
BEC  operators  with  Gaussian- shaped  correlation  functions 
and  their  approximations  (e.g.  Weaver  etai,  2003;  Di 
Lorenzo  et  aL,  2007).  Among  the  reasons  for  that  limited 
applicability  is  poor  conditioning  of  the  BEC  operators 
generated  by  high-degree  polynomials  and  the  necessity 
to  link  polynomial  coefficients  with  the  shape  of  the  BEC 
spectrum.  In  the  recent  studies  of  Hristopulos  and  Elogne 
(2007,  2009)  and  Yaremchuk  and  Smith  (2011),  correlation 
functions  associated  with  an  arbitrary  quadratic  polynomial 
of  the  homogeneous  diffusion  operator  were  obtained  and 
relationships  between  the  polynomial  coefficients  and  the 
magnitude/length  scale  of  the  corresponding  spectral  peak 
have  been  provided. 

In  this  note  the  result  of  Yaremchuk  and  Smith  (2011)  is 
extended  for  the  case  of  an  arbitrary  polynomial,  generating 
a  multiple-peak  BEC  spectrum.  Besides,  it  is  shown  that 
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the  action  of  the  BEC  operator  can  be  reduced  to  a 
sequence  of  inversions  of  the  quadratic  functions  of  the 
diffusion  operator,  thereby  relaxing  the  above-mentioned 
conditioning  problem. 

2.  Homogeneous  multi-scale  correlation  functions 

This  note  deals  with  covariance  modelling  in  Rw,  although 
the  results  can  be  extended  to  an  arbitrary  differentiable 
manifold  of  constant  curvature. 

A  general  form  of  the  inverse  BEC  operator  built  as  a 
polynomial  of  the  homogeneous  diffusion  operator  A  is 

L 

B"'  =I  +  £a,A'.  (I) 

/-I 

Here  I  is  the  identity  operator  and  ct{  are  real  numbers, 
constrained  by  the  positive  definiteness  requirement  of 
B“l,  which  can  be  taken  into  account  explicitly  by 
diagonalizing  B  via  the  Fourier  transform.  In  the 
Fourier  representation  the  inverse  BEC  operator  acts 
as  multiplication  by  the  polynomial  in  k2  *  |k|2  (k  is 
the  wavenumber),  and  the  positive-definiteness  property 
translates  into  the  requirement  for  the  spectral  polynomial 


and  can  be  interpreted  as  the  scales  (b~{)  and  magnitudes 
{a~] )  of  the  modes  forming  the  spectrum. 

Since  the  reciprocal  of  IT  '(Jt)  provides  the  spectral 
representation  of  the  BEC  operator,  the  matrix  elements 
B(r)  of  B  (covariance  functions)  depend  only  on  the  distance 
r  from  the  diagonal  and  can  be  obtained  in  the  form  of  a 
single  integral  over  k. 


B"(r)  = 


Zr~s  7  k*+lJs(kr)dk 
(2 ^jfJ  JWT^WTZ)- 

0  m 


(6) 


Here  /  denotes  the  Bessel  function  of  the  first  kind,  «  is  the 
dimension  of  the  physical  space  and  s  =  n/2  —  1.  Equation 
(6)  is  obtained  by  substitution  of  the  reciprocal  of  Eq. 
(3)  into  the  integral  of  the  inverse  Fourier  transform  and 
integrating  over  the  solid  angle  in  the  wavenumber  space 
(e.g.  Yaremchuk  and  Smith,  2011).  The  integral  (Eq.  (6)) 
can  be  taken  by  decomposing 


B(*)  = 


Z 

ri(*j+4K^+^) 


(7) 


into  elementary  fractions: 


L 

B-'(k2)  =  \  +  ']Tai(-k2)1  (2) 

l«l 


B(k)  = 


+ 


(8) 


to  be  positive  for  all  k2  >  0.  This  constraint  is  equivalent  to 
the  statement  that  the  right-hand  side  of  Eq.  (2)  must  not 
have  real  positive  roots.  A  particular  form  of  the  even-order 
polynomial  satisfying  this  requirement  is 

1  M 

-  rj(*i+^)(*2+£).  o) 

m  =  l 

where  M  =  E/2, 

^  =  <4> 

m 

where  overline  in  Eq.  (3)  denotes  complex  conjugate 
and  zm  —  am  4-  ibm  are  arbitrary  complex  numbers  with 
ambm  ^  0.  In  its  general  form,  the  polynomial  (3)  is 
additionally  multiplied  by  the  product  over  the  arbitrary 
number  of  real  negative  roots.  To  simplify  the  formulas,  we 
consider  this  case  in  the  Appendix,  and  focus  on  the  analysis 
of  Eq.  (3)  omitting  the  product  (summation)  limits  over 
m  and  assuming  there  are  no  real  negative  and  multiple 
roots.  The  latter  requirement  is  not  restrictive  for  practical 
purposes,  because  location  of  the  roots  is  never  known 
exactly,  and  the  BEC  spectrum  can  be  well  approximated  by 
Eq.  (3)  (see  Appendix). 

It  is  instructive  to  note  that  expression  (3)  can  also  be 
rewritten  in  the  form 

B_l  =  y  ril"">  +  ^  ~  +  (*  +  Bm)2|.  (5) 

m 

Compared  to  the  spectral  representation  (Eq.  (2)), 
representation  (5)  has  the  advantage  that  its  free  parameters 
are  not  constrained  by  the  positive-definiteness  requirement 


where 


—  .  - 


z 


(9) 


By  replacing  Eq.  (7)  in  Eq.  (6)  with  the  sum  (8)  the 
integral  is  reduced  to  the  sum  of  Hankel-Nicholson  type 
integrals  (e.g.  Abramowitz  and  Stegun,  1972,  eq.  11.4.44) 
and  can  be  taken  explicitly.  The  cited  integral  is  valid  if  the 
following  constraints  are  satisfied:  (a)  —1  <  $  <  3/2  and  (b) 
the  real  part  of  zm  is  positive.  The  first  condition  is  met  in  the 
practical  cases  of  0  <  m  <  5.  Furthermore,  in  view  of  Eq.  (5) 
both  am  and  bm  can  be  assumed  to  be  positive  without  loss 
of  generality  and  thus  the  second  constraint  is  also  satisfied. 
The  resulting  expression  for  #(r)  is 


B*(r)  = 


2r2~n 

(2*)? 


(tfmPrn  (Pm  ))  > 


(10) 


where  pm  =  zmrt  K  stands  for  the  modified  Bessel  function 
of  the  second  kind,  and  angle  brackets  denote  taking  the  real 
part. 

The  corresponding  correlation  functions  C*(r)  are 
obtained  through  normalizing  Eq.  ( 10)  by  B"( 0).  For  n  <  4 
the  BEC  function  values  at  r  =  0  are 


B‘(  0) 

—  y  .  (^m An) Urn i  2» 

rrt 

(ID 

b2(0) 

=  “  V'fam  logzm)> 

77  — 
m 

02) 

fl’(O) 

rrt 

(13) 

Copyright  CO  2012  Royal  Meteorological  Society 


Q.  f.  R  MrtroroL  Soc  138:  1948-1953  (2012) 


1950 


M.  Yaremchuk  and  A.  Scntchcv 


In  one-  and  three-dimensional  cases  the  correlation 
functions  can  also  be  expressed  In  terms  of  the  exponents: 


Clm  I 'A<lmZn,e  ^)/\Zm\2 

U^Z„)/\Z,„\2 

(14) 

<:V)=  1 

27T -r  ^2(qmzm) 

(15) 

for  Zm  to  initialize  an  iterative  procedure  of  approximating 
experimental  data. 

After  the  model  parameters  are  established,  the  action 
of  the  inverse  BEC  operator  can  be  computed  recursively 
by 

B“,  =  n[,-l^l'2A(2(4)I- A)].  (20) 

Ttt 


Relationships  ( 1 0)— ( 1 3)  provide  analytical  expressions 
for  the  multi-scale  homogeneous  BEC  functions  and  the 
corresponding  CFs.  In  many  applications,  it  is  often 
important  to  know  the  value  N  of  the  convolution  of  B  with 
the  <$ -function  at  r  =  0  (the  normalization  factor)  which  is 
used  for  constructing  the  BEC  operator  numerically.  The 
factor  can  be  found  by  integrating  C”(r)  over  R": 


2  y' 

B"(0)  ^  ‘ 


(lb) 


3.  Practical  issues 


In  applications,  a  BEC  model  is  often  constructed  by 
fitting  spectral  (7)  or  correlation  (10)  functions  to  those 
derived  from  experimental  data.  The  random  fields  under 
consideration  are  characterized  by  2 m  parameters,  which 
give  enough  freedom  for  approximating  complex  spectra. 
The  approximation  procedure  can  be  formulated  as  a  least 
squares  problem  in  2m  dimensions,  which  may  be  rather 
complicated  due  to  nonlinearity  of  B  with  respect  to  the 
fitting  parameters  am  and  bm .  Therefore  it  is  useful  to 
have  guidance  on  how  the  magnitudes  and  locations  of  the 
model  peaks  are  related  to  the  scales  and  amplitudes  of  the 
physical  modes  contributing  to  the  experimental  spectrum 
(Figure  1). 

The  contribution  of  the  mth  mode  to  the  spectrum  can 
be  assessed  by  integrating  the  right-hand  side  of  Eq.  (8): 


+ 


Hm 

*2  +  3,J 


j,  71 

dK  = - - — 

|Zm|2 


(17) 


In  the  limit  when  distances  \bi  —  bm\  between  the  spectral 
peaks  of  B  arc  much  larger  than  their  half-widths  am  (i.e. 
am/b„,  <K  0  in  particular),  Eq.  (17)  can  be  simplified  using 
the  asymptotic  approximations 


!  ib„ 


</m 


trl 


4wmn„ 


n^no-O*,2)2. 


so  that 


(18) 


Asymptotic  values  of  spectral  density  at  the  peaks  are 
respectively 


B(bm) 


tin 

4amn 


(19) 


i.e.  the  peak  amplitudes  are  inversely  proportional  to  a2m 
and  to  the  square  of  the  mode  scale  b~l.  Expressions 
(17M19)  can  be  useful  in  generating  the  first-guess  values 


The  inverse  BEC  model  (£q.  (20))  can  then  be  employed 
either  to  compute  the  action  of  B  using  an  iterative  inversion 
scheme  or  to  directly  compute  the  gradient  of  a  3dVar  cost 
function  involving  the  quadratic  form  xTB~lx,  where  x  is 
the  state  vector. 

The  considered  multi-scale  BEC  operators  can  be 
used  in  many  oceanographic  applications,  where  the 
background  errors  have  multi-scale  spectra.  For  instance, 
surface  waves  are  often  characterized  by  two-peak  spectra 
generated  by  the  swell  and  locally  forced  wind  waves. 
Filtering  such  a  wave-induced  signal  from  observations  ts 
important  in  many  applications  (e.g.  vertical  positioning 
of  the  autonomous  underwater  vehicles,  turbulence 
microstructure  measurements  in  shallow  seas). 

Figure  2(a)  demonstrates  typical  velocity  spectra, 
derived  from  observations  by  an  upward -looking  bottom - 
mounted  acoustic  Doppler  current  profiler  (Korotenko 
rtaL ,  2012).  Measurements  were  taken  in  the  period  of 
well -developed  wind  waves  with  a  dominant  frequency 
/  ~  0.2  Hz  superimposed  on  the  swell  (f  ^0.1  Hz) 
propagating  from  the  Bay  of  Biscay.  Slight  asymmetry 
of  the  beam  directions  with  respect  to  the  vertical 
prevents  cancellation  of  the  wave-induced  orbital  motions 
in  averaging  over  the  beams  and  contaminates  the 
turbulence  spectrum  with  a  double-bump  feature  seen  in 
Figure  2(a). 

Impact  of  the  surface  waves  can  be  removed  by 
constructing  a  rational  filter  (e.g  Antoniou,  2000)  F(k)  = 
H(A)B“ 1  (A)  using  the  polynomials  (1)  of  the  diffusion 
operator  A  =  d^.  The  rational  function  F(k)  can  be  obtained 
by  adjusting  the  filter  parameters  z,  z  to  the  ratio  between 
the  observed  spectrum  and  its  power-law  approximation 
(dashed  line  in  Figure  2(a))  in  the  wave-contaminated 
frequency  band  co  =  [0.06  -  0.4]  Hz.  After  the  adjustment, 
the  matrix  elements  of  B  are  computed  using  Eqs  (10), 
(14): 


m-l 


7^7  cos|frmf  +  arg(qmzm)|, 


(21) 


whereas  the  action  of  B~ 1  is  given  by  Eq. 

(20): 


=  [I  I1  -  l^l'23«(2(^)  -  *>«>]  •  (22) 

The  filter  is  implemented  by  differentiating  the  series 
with  Eq,  (22)  and  then  smoothing  it  with  the  kernel 
(Eq.  (21 )). 

Figure  2(b)  demonstrates  the  result  of  fitting  the  filter 
parameters  zmt  zm  to  the  wave- induced  part  of  the  spectrum 
shown  in  Figure  2(a),  The  fit  has  a  relative  error  of  7%  within 
the  target  band  of  0.06-0.4  Hz.  In  the  same  frequency  band, 
the  filtered  series  spectrum  has  similar  (6%)  deviation  from 
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Figure  I.  An  example  oflhc  normalized  spectrum  B(Jt)  (Eq.  (7),  left  panel)  and  the  respective  correlation  function  fl(r)//?(0)  (Eqs  (10>  12).  right  panel) 
for  M  *=  n  =  2,  Z|  «  .5  *f  3i;  zj  =  .2  4-  6i. 


by  surface  waves  is  seen  as  a  pronounced  double  peak.  Grey  line  shows  the  spectrum  of  filtered  series,  (b)  Solid  black  line  b  the  ratio  between  the  observed 
spectrum  and  its  power  law  interpolation  (dashed  line  in  left  pand)  within  the  wave-contaminated  band  a>  This  ratio  was  used  for  constructing  the  filler 
¥,  whose  spectrum  is  shown  in  grey 


the  power  law  when  normalized  by  deviation  of  the  observed 
spectrum  from  that  law.  Proximity  of  the  filtered  spectrum 
to  the  power  law  can  be  further  improved  by  increasing  the 
order  M  of  the  polynomials  B~'(k)  and  B~l( k). 

4.  Summary  and  discussion 

Analytical  expressions  for  the  matrix  elements  of  homo¬ 
geneous  BEC  operators  generated  by  the  polynomials 
(Eq.  (1))  of  the  diffusion  operator  are  obtained.  The 
considered  BEC  operators  can  be  used  in  geophysical 
applications  involving  multi- scale  phenomena  whose  con¬ 
tribution  to  the  spectrum  can  be  modelled  by  adjusting 
the  free  parameters  (polynomial  coefficients)  of  the  BEC 
model.  Applicability  of  the  technique  to  a  simple  two- 
scale  one-dimensional  filtering  problem  has  been  demon¬ 
strated. 

A  particular  advantage  of  the  considered  type  of  BEC 
operators  is  the  fact  that  their  inverses  can  be  represented 
by  sparse  matrices  that  can  be  efficiently  implemented 
on  the  grids  of  various  complexity.  Explicit  partitioning 
of  the  inverse  operators  (Eq.  (3))  ensures  their  positive- 
definiteness  and  provides  a  recursive  algorithm  (Eq.  (20)) 
for  computing  the  action  of  the  BEC  operator  which  has 
reasonably  conditioned  matrices  on  each  iteration. 
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Presented  results  arc  also  valid  in  the  homogeneous 
anisotropic  case,  because  the  latter  can  be  reduced  to 
isotropic  form  by  the  appropriate  coordinate  transformation 
(e.g.  Xu,  2005;  Yaremchuk  and  Carrier,  2012).  The 
obtained  analytical  expressions  for  the  correlation  functions 
(Eqs  (10-15))  can  be  useful  in  finding  the  BEC  model 
parameters  for  the  fields  whose  local  decorrelation  scales 
p  do  not  change  significantly  at  distances  of  the  order 
of  p. 

In  the  more  general  inhomogeneous  case  analytical 
formulas  for  the  matrix  elements  of  B  cannot  be  obtained, 
and  inversion  of  the  operator  (Eq.  (20))  has  to  be  performed 
numerically.  For  the  BEC  models  with  M  >  2  such  inversion 
may  encounter  difficulties  associated  with  the  condition 
number  of  B,  which  grows  exponentially  with  the  number  of 
model  parameters.  In  view  of  the  decomposition  (Eq.  (20)), 
however,  this  inversion  can  be  performed  consecutively 
by  iterative  solutions  of  M  linear  systems  whose  condition 
numbers  are  limited  from  above  by  the  maximum  eigenvalue 
<>fi£r2A2. 

In  higher  dimensions  (n  >  1)  the  polynomial  BEC 
model  can  be  further  improved  by  introducing  anisotropic 
inhomogeneous  diffusion  operators  separately  for  each 
mode.  The  respective  diffusion  tensors  can  be  adjusted  using 
prior  knowledge  of  the  impact  of  the  background  fields  on 
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Figure  3.  An  example  of  the  normal ized  spectrum  generated  by  adding 
two  negative  roots  p\  =  I  and  fo  =  4  to  the  spectrum  in  Figure  1 
Ui  =  .5  +  3i:  Zi  =  .2  +  6i).  Grey  spect rum  is  the  approximation  obtained 
by  replacing  pi  and  pj  with  z\  — 1.46  +  .01  i.  The  approximation  error  is 

5%. 


the  error  characteristics  described  by  the  corresponding 
spectral  peak. 

For  the  BEC  spectra  characterized  by  deep  gaps  between 
the  peaks,  the  multi-scale  approach  of  Li  etal  (pers. 
comm.  2012)  may  prove  to  be  more  computationally 
effective,  as  it  does  not  take  into  the  account  scale 
interactions,  and  adopts  the  additive  BEC  model.  In  the 
case  when  a  pair  of  closely  spaced  peaks  (Figure  2(a)) 
exists,  the  technique  of  Li  etal  (pers.  comm.  2012)  can 
be  generalized  by  introducing  a  two-parameter  model 
to  account  for  the  additive  BEC  component  at  the 
corresponding  scale. 

Appendix 


Real  negative  roots  provide  limited  freedom  to  controlling 
l  he  shape  of  the  BEC  spectrum  because  in  this  case  the  poles 
of  Eq.  (7)  are  located  outside  the  range  k  >  0.  Therefore, 
these  poles  just  add  a  monotonously  decaying  function  of  k 
which  has  the  largest  impact  on  the  long-wave  part  of  B(k). 
Furthermore,  spectral  contribution  of  the  negative  roots  can 
always  be  well  approximated  by  a  pair  of  complex  roots 
e2  —  b2  ±  libs,  where  s  is  a  small  number  (Figure  3). 

Taking  N  =  L  —  2 M  negative  roots  -p„  (pH  >  0)  into 
account  supplements  all  the  formulas  with  an  extra 
summation  (product)  over  these  roots.  In  the  following, 
the  ‘negative-root  generalizations'  of  the  key  formulas  are 
listed.  For  clarity,  we  keep  the  numbering  and  abbreviate 
sums/products  from  the  main  text  by  (E)  and  |F1} 
respectively: 

= ^n*ri(Ar2 <a3> 

l 

*-<n»n*-<nin  Pn.  <A4) 

«=2M+1  n 


(jm  — 


‘b,  = 


z 


(*m-*m>(n!n(*m+P«) 

ti 

z 


;  0  <  m  <  M, 


n  ip«+4itkp»-p/) 

m  j+n 


;  1M  <n<L,  (A8) 


B"(r)= 


1 

(2jt)5 


2<Ei  +  'Ew,mp*k' 

n 


where  p„  =  p„r,  and  p.  =  Jp~n : 
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